Practical 4: Multispecies
Olivier Gimenez
May 2023
- 1 Dapprima selezione delle variabili per ciascuna delle quattro specie
- 2 Data
- 3 Single species
- 4 Test delle ipotesi ecologiche di
interazione inter-specifica
- 4.1 Interazione stambecco - bestiame domestico
- 4.2 Interazione leopardo - bestiame domestico
- 4.3 Interazione lupo - bestiame domestico
- 4.4 Interazione leopardo - stambecco
- 4.5 Interazione leopardo - lupo
- 4.6 Interazione stambecco - lupo
- 4.7 Riassunto grafico della selezione dei modelli per le interazioni inter-specifiche
1 Dapprima selezione delle variabili per ciascuna delle quattro specie
2 Data
Thanks Marco Salvatori.
covar<- read.csv2('data/covs_original_ALL.csv', stringsAsFactors = F)
dt1<- covar[,c('Sampling.Unit', 'Area', 'CT.sens', 'EL', 'SL', 'Dis')]
dt1[,4:6] <- apply(X=dt1[,4:6], MARGIN=2, FUN=decostand, method='standardize')
var<- dt1ll<-readRDS(file = 'data/matrici_complessive.rds')dom1<-ll$domestici
dom<-ll$domestici
ibex<-ll$ibex
leo<-ll$leo
lupus<-ll$lupo3 Single species
DOMESTICI
do<- unmarkedFrameOccu(y=dom, siteCovs = var, obsCovs = NULL)
ib<- unmarkedFrameOccu(y=ibex, siteCovs = var, obsCovs = NULL)
le<- unmarkedFrameOccu(y=leo, siteCovs = var, obsCovs = NULL)
lu<- unmarkedFrameOccu(y=lupus, siteCovs = var, obsCovs = NULL)3.1 Bestiame
##p
mod0<-occu(~1~Area+EL+SL+Dis ,do)
mod1<-occu(~Dis~Area+EL+SL+Dis,do)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,do)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, do)mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)
modSel(mods)## Hessian is singular.## nPars AIC delta AICwt cumltvWt
## Dis+CT.sens 10 5230.07 0.00 1.0e+00 1.00
## CT.sens 9 5317.00 86.93 1.3e-19 1.00
## Zero 8 5363.11 133.03 1.3e-29 1.00
## Dis 9 5832.17 602.09 1.8e-131 1.00## psi
mod1<-occu(~Dis+CT.sens~1, do)
mod2<-occu(~Dis+CT.sens~Area, do)
mod3<-occu(~Dis+CT.sens~Dis, do)
mod4<-occu(~Dis+CT.sens~EL, do)
mod5<-occu(~Dis+CT.sens~SL, do)
mod6<-occu(~Dis+CT.sens~Area+SL, do)
mod7<-occu(~Dis+CT.sens~Area+EL, do)
mod8<-occu(~Dis+CT.sens~Area+Dis, do)
mod9<-occu(~Dis+CT.sens~EL+Dis, do)
mod10<-occu(~Dis+CT.sens~EL+SL, do)
mod11<-occu(~Dis+CT.sens~Dis+SL, do)
mod12<-occu(~Dis+CT.sens~Area+EL+SL, do)
mod13<-occu(~Dis+CT.sens~Area+EL+Dis, do)
mod14<-occu(~Dis+CT.sens~Area+SL+Dis, do)
mod15<-occu(~Dis+CT.sens~SL+EL+Dis, do)
mod16<-occu(~Dis+CT.sens~ Area+EL+SL+Dis, do)mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7,
'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13, 'Area+SL+Dis'=mod14,
'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)ms<-modSel(mdls)msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms## model negLogLike nPars AIC delta AICwt
## 14 Area+EL+Dis 2605.266 9 5228.533 0.000000 6.755202e-01
## 17 Area+EL+SL+Dis 2605.036 10 5230.072 1.539166 3.129051e-01
## 8 Area+EL 2611.004 8 5238.007 9.474579 5.919155e-03
## 13 Area+EL+SL 2610.830 9 5239.660 11.127025 2.590810e-03
## 10 EL+Dis 2614.273 6 5240.545 12.012452 1.664055e-03
## 16 SL+EL+Dis 2614.034 7 5242.069 13.536115 7.767984e-04
## 5 EL 2616.613 5 5243.226 14.693349 4.355301e-04
## 11 EL+SL 2616.469 6 5244.939 16.406348 1.849461e-04
## 3 Area 2619.938 7 5253.875 25.342434 2.121283e-06
## 7 Area+SL 2619.828 8 5255.657 27.123998 8.704361e-07
## 12 Dis+SL 2622.515 6 5257.029 28.496660 4.381948e-07
## 2 K 2627.666 4 5263.332 34.799490 1.875103e-08
## 6 SL 2627.600 5 5265.199 36.666602 7.372025e-09
## 1 Zero 2673.553 8 5363.106 134.572937 4.050386e-30
## 15 Area+SL+Dis 2868.362 9 5754.725 526.192226 3.702100e-115
## 4 Dis 2904.134 5 5818.269 589.736298 5.888821e-129
## 9 Area+Dis 2904.189 8 5824.379 595.846176 2.775142e-130summary(mod13)##
## Call:
## occu(formula = ~Dis + CT.sens ~ Area + EL + Dis, data = do)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.131 0.277 0.471 0.63789
## AreaSB -0.688 0.416 -1.654 0.09803
## AreaSU 1.008 0.475 2.121 0.03389
## AreaTB -0.811 0.459 -1.765 0.07758
## EL -0.582 0.187 -3.111 0.00187
## Dis -0.617 0.190 -3.241 0.00119
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -2.436 0.0504 -48.34 0.00e+00
## Dis -0.467 0.0552 -8.47 2.53e-17
## CT.sensmedium 0.487 0.0850 5.73 1.02e-08
##
## AIC: 5228.533
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 87
## Bootstrap iterations: 03.2 Stambecco Siberiano
### p
mod0<-occu(~1~Area+EL+SL+Dis ,ib)
mod1<-occu(~Dis~Area+EL+SL+Dis,ib)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,ib)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, ib)mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)
modSel(mods)## nPars AIC delta AICwt cumltvWt
## Dis+CT.sens 10 1444.34 0.00 9.6e-01 0.96
## Dis 9 1450.47 6.13 4.5e-02 1.00
## Zero 8 1558.03 113.68 2.0e-25 1.00
## CT.sens 9 1560.50 116.16 5.7e-26 1.00## psi
mod1<-occu(~CT.sens~1, ib)
mod2<-occu(~CT.sens~Area, ib)
mod3<-occu(~CT.sens~Dis, ib)
mod4<-occu(~CT.sens~EL, ib)
mod5<-occu(~CT.sens~SL, ib)
mod6<-occu(~CT.sens~Area+SL, ib)
mod7<-occu(~CT.sens~Area+EL, ib)
mod8<-occu(~CT.sens~Area+Dis, ib)
mod9<-occu(~CT.sens~EL+Dis, ib)
mod10<-occu(~CT.sens~EL+SL, ib)
mod11<-occu(~CT.sens~Dis+SL, ib)
mod12<-occu(~CT.sens~Area+EL+SL, ib)
mod13<-occu(~CT.sens~Area+EL+Dis, ib)
mod14<-occu(~CT.sens~Area+SL+Dis, ib)
mod15<-occu(~CT.sens~SL+EL+Dis, ib)
mod16<-occu(~CT.sens~ Area+EL+SL+Dis, ib)mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7,
'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13, 'Area+SL+Dis'=mod14,
'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)ms<-modSel(mdls)## Hessian is singular.
## Hessian is singular.msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms## model negLogLike nPars AIC delta AICwt
## 13 Area+EL+SL 714.4724 8 1444.945 0.000000 4.711932e-01
## 6 SL 719.1642 4 1446.328 1.383615 2.359124e-01
## 11 EL+SL 718.6172 5 1447.234 2.289506 1.499820e-01
## 16 SL+EL+Dis 718.3004 6 1448.601 3.656000 7.573712e-02
## 5 EL 721.3001 4 1450.600 5.655432 2.787004e-02
## 4 Dis 721.3584 4 1450.717 5.771906 2.629333e-02
## 10 EL+Dis 721.0618 5 1452.124 7.178808 1.301190e-02
## 8 Area+EL 751.9379 7 1517.876 72.930922 6.861998e-17
## 3 Area 753.7809 6 1519.562 74.617022 2.953375e-17
## 14 Area+EL+Dis 751.9234 8 1519.847 74.901928 2.561251e-17
## 9 Area+Dis 753.6673 7 1521.335 76.389846 1.217180e-17
## 1 Zero 771.0131 8 1558.026 113.081478 1.311801e-25
## 17 Area+EL+SL+Dis 771.2512 9 1560.502 115.557665 3.803390e-26
## 15 Area+SL+Dis 778.7924 8 1573.585 128.640099 5.487271e-29
## 7 Area+SL 781.5926 7 1577.185 132.240404 9.069018e-30
## 2 K 794.2281 3 1594.456 149.511353 1.611467e-33
## 12 Dis+SL 794.2286 5 1598.457 153.512482 2.179652e-34summary(mod12)##
## Call:
## occu(formula = ~CT.sens ~ Area + EL + SL, data = ib)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.682 0.307 -2.221 0.0263
## AreaSB -0.576 0.490 -1.175 0.2399
## AreaSU 0.174 0.424 0.411 0.6812
## AreaTB -1.295 0.615 -2.106 0.0352
## EL 0.270 0.187 1.448 0.1476
## SL 0.368 0.169 2.173 0.0298
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -3.520 0.112 -31.55 1.92e-218
## CT.sensmedium 0.684 0.227 3.02 2.55e-03
##
## AIC: 1444.945
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 130
## Bootstrap iterations: 03.3 Leopardo delle nevi
### p
mod0<-occu(~1~Area+EL+SL+Dis ,le)
mod1<-occu(~Dis~Area+EL+SL+Dis,le)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,le)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, le)mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)#, 'Dis+CT.sens+Area'=mod4 )
modSel(mods)## nPars AIC delta AICwt cumltvWt
## Dis+CT.sens 10 2481.82 0.00 5.5e-01 0.55
## CT.sens 9 2482.22 0.41 4.5e-01 1.00
## Dis 9 2524.36 42.54 3.2e-10 1.00
## Zero 8 2531.31 49.50 9.8e-12 1.00## psi
mod1<-occu(~CT.sens~1, le)
mod2<-occu(~CT.sens~Area, le, starts = rep(-1,6))
mod3<-occu(~CT.sens~Dis, le)
mod4<-occu(~CT.sens~EL, le)
mod5<-occu(~CT.sens~SL, le)
mod6<-occu(~CT.sens~Area+SL, le, starts = rep(-1,7))
mod7<-occu(~CT.sens~Area+EL, le)
mod8<-occu(~CT.sens~Area+Dis, le, starts = rep(1,7))
mod9<-occu(~CT.sens~EL+Dis, le)
mod10<-occu(~CT.sens~EL+SL, le)
mod11<-occu(~CT.sens~Dis+SL, le)
mod12<-occu(~CT.sens~Area+EL+SL, le)
mod13<-occu(~CT.sens~Area+EL+Dis, le, starts = rep(1,8))
mod14<-occu(~CT.sens~Area+SL+Dis, le, starts = rep(1,8))
mod15<-occu(~CT.sens~SL+EL+Dis, le)
mod16<-occu(~CT.sens~ Area+EL+SL+Dis, le)mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7,
'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13, 'Area+SL+Dis'=mod14,
'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)ms<-modSel(mdls)msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms## model negLogLike nPars AIC delta AICwt
## 3 Area 1233.772 6 2479.543 0.00000000 2.089589e-01
## 9 Area+Dis 1232.798 7 2479.596 0.05234249 2.035611e-01
## 7 Area+SL 1232.951 7 2479.903 0.35945726 1.745845e-01
## 15 Area+SL+Dis 1232.205 8 2480.410 0.86671176 1.354742e-01
## 14 Area+EL+Dis 1232.715 8 2481.429 1.88585873 8.138639e-02
## 8 Area+EL 1233.769 7 2481.539 1.99546759 7.704607e-02
## 13 Area+EL+SL 1232.951 8 2481.903 2.35929613 6.423122e-02
## 17 Area+EL+SL+Dis 1232.111 9 2482.222 2.67844085 5.475764e-02
## 10 EL+Dis 1257.842 5 2525.685 46.14126918 1.998074e-11
## 16 SL+EL+Dis 1257.513 6 2527.026 47.48274234 1.021679e-11
## 1 Zero 1257.656 8 2531.312 51.76872549 1.198464e-12
## 4 Dis 1262.191 4 2532.382 52.83916398 7.017495e-13
## 12 Dis+SL 1261.839 5 2533.678 54.13488590 3.671298e-13
## 5 EL 1267.189 4 2542.378 62.83454751 4.739278e-15
## 2 K 1268.284 3 2542.569 63.02572991 4.307224e-15
## 11 EL+SL 1266.984 5 2543.968 64.42488156 2.139811e-15
## 6 SL 1268.066 4 2544.132 64.58845437 1.971769e-15summary(mod2)##
## Call:
## occu(formula = ~CT.sens ~ Area, data = le, starts = rep(-1, 6))
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.753 0.305 -2.47 1.36e-02
## AreaSB 1.805 1.599 1.13 2.59e-01
## AreaSU 2.293 0.465 4.93 8.36e-07
## AreaTB -2.139 0.789 -2.71 6.69e-03
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -3.05 0.0661 -46.16 0.000000
## CT.sensmedium -1.57 0.4430 -3.54 0.000405
##
## AIC: 2479.543
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 26
## Bootstrap iterations: 03.4 Lupo
### p
mod0<-occu(~1~Area+EL+SL+Dis ,lu)
mod1<-occu(~Dis~Area+EL+SL+Dis,lu)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,lu)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, lu)mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)#, 'Dis+CT.sens+Area'=mod4 )
modSel(mods)## nPars AIC delta AICwt cumltvWt
## Zero 8 1043.98 0.00 0.523 0.52
## CT.sens 9 1045.83 1.85 0.208 0.73
## Dis 9 1045.98 2.00 0.193 0.92
## Dis+CT.sens 10 1047.82 3.84 0.077 1.00## psi
mod1<-occu(~1~1, lu)
mod2<-occu(~1~Area, lu)
mod3<-occu(~1~Dis, lu)
mod4<-occu(~1~EL, lu)
mod5<-occu(~1~SL, lu)
mod6<-occu(~1~Area+SL, lu)
mod7<-occu(~1~Area+EL, lu)
mod8<-occu(~1~Area+Dis, lu)
mod9<-occu(~1~EL+Dis, lu)
mod10<-occu(~1~EL+SL, lu)
mod11<-occu(~1~Dis+SL, lu)
mod12<-occu(~1~Area+EL+SL, lu)
mod13<-occu(~1~Area+EL+Dis, lu)
mod14<-occu(~1~Area+SL+Dis, lu)
mod15<-occu(~1~SL+EL+Dis, lu)
mod16<-occu(~1~ Area+EL+SL+Dis, lu)mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7,
'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13, 'Area+SL+Dis'=mod14,
'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)ms<-modSel(mdls)msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms## model negLogLike nPars AIC delta AICwt
## 10 EL+Dis 515.8354 4 1039.671 0.000000 0.254151574
## 4 Dis 516.9341 3 1039.868 0.197321 0.230274105
## 16 SL+EL+Dis 515.7615 5 1041.523 1.852063 0.100675192
## 12 Dis+SL 516.8750 4 1041.750 2.079139 0.089869740
## 14 Area+EL+Dis 514.1072 7 1042.214 2.543562 0.071246811
## 2 K 519.6448 2 1043.290 3.618671 0.041620605
## 9 Area+Dis 515.7373 6 1043.475 3.803640 0.037943976
## 1 Zero 513.9901 8 1043.980 4.309232 0.029468296
## 17 Area+EL+SL+Dis 513.9901 8 1043.980 4.309232 0.029468296
## 3 Area 517.1861 5 1044.372 4.701246 0.024223126
## 5 EL 519.3935 3 1044.787 5.116011 0.019686343
## 8 Area+EL 516.6010 6 1045.202 5.531086 0.015996782
## 15 Area+SL+Dis 515.6164 7 1045.233 5.561959 0.015751745
## 6 SL 519.6207 3 1045.241 5.570524 0.015684436
## 7 Area+SL 517.0778 6 1046.156 6.484662 0.009930393
## 11 EL+SL 519.3674 4 1046.735 7.064001 0.007433007
## 13 Area+EL+SL 516.4900 7 1046.980 7.309139 0.006575572summary(mod9)##
## Call:
## occu(formula = ~1 ~ EL + Dis, data = lu)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.764 0.221 -3.46 0.000535
## EL -0.299 0.205 -1.46 0.143235
## Dis 0.506 0.198 2.55 0.010718
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -3.95 0.143 -27.6 8.16e-168
##
## AIC: 1039.671
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 49
## Bootstrap iterations: 03.5 Synthesis
Perciò i migliori modelli di singola specie, selezionando il migliore in termini di AIC, sono i seguenti:
## p psi
## leopardo ~CT.sens ~Area
## lupo ~1 ~Elevation+Distance
## stambecco ~CT.sens ~Area+Elevation+Slope
## bestiame Distance+CT.sens ~Area+Elevation+Distance3.6 Occupancy di singola specie - grafici
4 Test delle ipotesi ecologiche di interazione inter-specifica
4.1 Interazione stambecco - bestiame domestico
dom_ibex_dip <- occuMulti(detformulas = c('~Dis+CT.sens','~CT.sens'), stateformulas = c('~Area+EL+Dis', '~Area+EL+SL', '~1'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))
dom_ibex_ind <- occuMulti(detformulas = c('~Dis+CT.sens','~CT.sens'), stateformulas = c('~Area+EL+Dis', '~Area+EL+SL', '~0'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))library(AICcmodavg)
modavgShrink(list(dom_ibex_dip, dom_ibex_ind), parm='[domestici:ibex] (Intercept)', parm.type='psi')## model negLogLike nPars AIC delta AICwt
## 2 dom_ibex_ind 3319.739 17 6673.477 0.00000000 0.5093794
## 1 dom_ibex_dip 3318.776 18 6673.552 0.07504392 0.4906206Piccola differenza, ma il modello con l’AIC migliore include l’interazione domestici-stambecco.
Modello senza interazione:
##
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis",
## "~Area+EL+SL", "~0"), data = dt1, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) 0.130 0.277 0.471 0.63796
## [domestici] AreaSB -0.687 0.416 -1.654 0.09816
## [domestici] AreaSU 1.008 0.475 2.121 0.03389
## [domestici] AreaTB -0.811 0.459 -1.766 0.07748
## [domestici] EL -0.582 0.187 -3.111 0.00187
## [domestici] Dis -0.617 0.190 -3.241 0.00119
## [ibex] (Intercept) -0.681 0.307 -2.220 0.02639
## [ibex] AreaSB -0.576 0.490 -1.175 0.23987
## [ibex] AreaSU 0.174 0.424 0.411 0.68141
## [ibex] AreaTB -1.295 0.615 -2.106 0.03517
## [ibex] EL 0.270 0.187 1.448 0.14756
## [ibex] SL 0.368 0.169 2.173 0.02980
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) -2.436 0.0504 -48.34 0.00e+00
## [domestici] Dis -0.467 0.0552 -8.47 2.53e-17
## [domestici] CT.sensmedium 0.487 0.0850 5.73 1.02e-08
## [ibex] (Intercept) -3.520 0.1116 -31.55 1.97e-218
## [ibex] CT.sensmedium 0.684 0.2266 3.02 2.54e-03
##
## AIC: 6673.477
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 88
## Bootstrap iterations: 0Modello con interazione:
##
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis",
## "~Area+EL+SL", "~1"), data = dt1, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) 0.318 0.312 1.020 0.30796
## [domestici] AreaSB -0.766 0.423 -1.809 0.07049
## [domestici] AreaSU 1.027 0.479 2.146 0.03190
## [domestici] AreaTB -0.938 0.472 -1.988 0.04682
## [domestici] EL -0.564 0.189 -2.992 0.00277
## [domestici] Dis -0.616 0.190 -3.239 0.00120
## [ibex] (Intercept) -0.398 0.370 -1.076 0.28198
## [ibex] AreaSB -0.687 0.501 -1.372 0.16994
## [ibex] AreaSU 0.234 0.430 0.544 0.58620
## [ibex] AreaTB -1.377 0.620 -2.221 0.02632
## [ibex] EL 0.184 0.198 0.930 0.35230
## [ibex] SL 0.366 0.169 2.169 0.03006
## [domestici:ibex] (Intercept) -0.522 0.379 -1.376 0.16867
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) -2.437 0.0504 -48.32 0.00e+00
## [domestici] Dis -0.467 0.0552 -8.47 2.53e-17
## [domestici] CT.sensmedium 0.487 0.0850 5.73 1.01e-08
## [ibex] (Intercept) -3.522 0.1118 -31.51 6.95e-218
## [ibex] CT.sensmedium 0.687 0.2263 3.04 2.39e-03
##
## AIC: 6673.552
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 152
## Bootstrap iterations: 04.2 Interazione leopardo - bestiame domestico
## model negLogLike nPars AIC delta AICwt
## 1 dom_leo_dip 3837.190 16 7706.380 0.000000 0.7496418
## 2 dom_leo_ind 3839.286 15 7708.573 2.193406 0.2503582Il modello migliore per AIC contiene l’interazione domestici-leopardo.
Modello senza interazione:
##
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis",
## "~Area", "~0"), data = dt2, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) 0.130 0.277 0.469 6.39e-01
## [domestici] AreaSB -0.689 0.416 -1.657 9.74e-02
## [domestici] AreaSU 1.009 0.475 2.124 3.37e-02
## [domestici] AreaTB -0.810 0.459 -1.763 7.79e-02
## [domestici] EL -0.582 0.187 -3.112 1.86e-03
## [domestici] Dis -0.616 0.190 -3.239 1.20e-03
## [leopardo] (Intercept) -0.724 0.303 -2.388 1.70e-02
## [leopardo] AreaSB 9.017 NaN NaN NaN
## [leopardo] AreaSU 2.263 0.464 4.877 1.08e-06
## [leopardo] AreaTB -2.167 0.789 -2.746 6.04e-03
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) -2.436 0.0504 -48.34 0.00e+00
## [domestici] Dis -0.467 0.0552 -8.47 2.51e-17
## [domestici] CT.sensmedium 0.487 0.0850 5.73 1.03e-08
## [leopardo] (Intercept) -3.052 0.0661 -46.15 0.00e+00
## [leopardo] CT.sensmedium -1.840 0.2427 -7.58 3.39e-14
##
## AIC: 7708.573
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 75
## Bootstrap iterations: 0Modello con interazione:
##
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis",
## "~Area", "~1"), data = dt2, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) 0.4706 0.340 1.383 1.67e-01
## [domestici] AreaSB -0.0494 0.588 -0.084 9.33e-01
## [domestici] AreaSU 1.4974 0.548 2.731 6.31e-03
## [domestici] AreaTB -1.1095 0.495 -2.242 2.50e-02
## [domestici] EL -0.5819 0.187 -3.110 1.87e-03
## [domestici] Dis -0.6140 0.190 -3.232 1.23e-03
## [leopardo] (Intercept) -0.1707 0.418 -0.408 6.83e-01
## [leopardo] AreaSB 6.2365 54.641 0.114 9.09e-01
## [leopardo] AreaSU 2.4152 0.491 4.919 8.70e-07
## [leopardo] AreaTB -2.5134 0.821 -3.061 2.21e-03
## [domestici:leopardo] (Intercept) -0.9819 0.507 -1.935 5.30e-02
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) -2.437 0.0505 -48.30 0.00e+00
## [domestici] Dis -0.467 0.0552 -8.47 2.51e-17
## [domestici] CT.sensmedium 0.487 0.0850 5.73 1.02e-08
## [leopardo] (Intercept) -3.053 0.0661 -46.16 0.00e+00
## [leopardo] CT.sensmedium -1.834 0.3154 -5.81 6.06e-09
##
## AIC: 7706.38
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 85
## Bootstrap iterations: 04.3 Interazione lupo - bestiame domestico
## model negLogLike nPars AIC delta AICwt
## 1 dom_lup_dip 3118.852 14 6265.704 0.000000 0.7772587
## 2 dom_lup_ind 3121.102 13 6268.203 2.499525 0.2227413Anche in questo caso il modello con migliore AIC prevede l’interazione.
Modello senza interazione:
##
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~1"), stateformulas = c("~Area+EL+Dis",
## "~EL+Dis", "~0"), data = dt3, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) 0.131 0.277 0.471 0.637865
## [domestici] AreaSB -0.688 0.416 -1.655 0.098020
## [domestici] AreaSU 1.008 0.475 2.121 0.033895
## [domestici] AreaTB -0.811 0.459 -1.765 0.077574
## [domestici] EL -0.582 0.187 -3.111 0.001867
## [domestici] Dis -0.617 0.190 -3.241 0.001190
## [lupo] (Intercept) -0.764 0.221 -3.462 0.000537
## [lupo] EL -0.299 0.205 -1.463 0.143350
## [lupo] Dis 0.506 0.198 2.552 0.010715
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) -2.436 0.0504 -48.34 0.00e+00
## [domestici] Dis -0.467 0.0552 -8.47 2.53e-17
## [domestici] CT.sensmedium 0.487 0.0850 5.73 1.02e-08
## [lupo] (Intercept) -3.950 0.1431 -27.61 8.46e-168
##
## AIC: 6268.203
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 118
## Bootstrap iterations: 0Modello con interazione:
##
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~1"), stateformulas = c("~Area+EL+Dis",
## "~EL+Dis", "~1"), data = dt3, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) -0.154 0.310 -0.495 0.620390
## [domestici] AreaSB -0.688 0.416 -1.656 0.097648
## [domestici] AreaSU 1.002 0.475 2.111 0.034794
## [domestici] AreaTB -0.813 0.460 -1.768 0.077043
## [domestici] EL -0.548 0.191 -2.873 0.004071
## [domestici] Dis -0.733 0.205 -3.579 0.000345
## [lupo] (Intercept) -1.300 0.341 -3.807 0.000140
## [lupo] EL -0.187 0.209 -0.893 0.372037
## [lupo] Dis 0.599 0.210 2.848 0.004406
## [domestici:lupo] (Intercept) 0.912 0.441 2.068 0.038680
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [domestici] (Intercept) -2.436 0.0504 -48.35 0.00e+00
## [domestici] Dis -0.467 0.0552 -8.47 2.55e-17
## [domestici] CT.sensmedium 0.487 0.0850 5.72 1.05e-08
## [lupo] (Intercept) -3.939 0.1418 -27.77 1.06e-169
##
## AIC: 6265.704
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 144
## Bootstrap iterations: 04.4 Interazione leopardo - stambecco
## model negLogLike nPars AIC delta AICwt
## 1 ib_leo_dip 1944.109 15 3918.218 0.000000 0.9583079
## 2 ib_leo_ind 1948.244 14 3924.488 6.269715 0.0416921In questo caso il modello con AIC più basso prevede l’interazione. Modello senza interazione:
##
## Call:
## occuMulti(detformulas = c("~CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+SL",
## "~Area", "~0"), data = dt4, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [ibex] (Intercept) -0.682 0.307 -2.221 2.63e-02
## [ibex] AreaSB -0.576 0.490 -1.175 2.40e-01
## [ibex] AreaSU 0.174 0.424 0.411 6.81e-01
## [ibex] AreaTB -1.295 0.615 -2.106 3.52e-02
## [ibex] EL 0.270 0.187 1.448 1.48e-01
## [ibex] SL 0.368 0.169 2.172 2.98e-02
## [leopardo] (Intercept) -0.753 0.305 -2.469 1.36e-02
## [leopardo] AreaSB 1.790 1.569 1.140 2.54e-01
## [leopardo] AreaSU 2.291 0.465 4.926 8.40e-07
## [leopardo] AreaTB -2.142 0.789 -2.713 6.66e-03
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [ibex] (Intercept) -3.520 0.1116 -31.55 1.93e-218
## [ibex] CT.sensmedium 0.684 0.2266 3.02 2.54e-03
## [leopardo] (Intercept) -3.052 0.0661 -46.17 0.00e+00
## [leopardo] CT.sensmedium -1.565 0.4402 -3.55 3.79e-04
##
## AIC: 3924.488
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 109
## Bootstrap iterations: 0Modello con interazione:
##
## Call:
## occuMulti(detformulas = c("~CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+SL",
## "~Area", "~1"), data = dt4, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [ibex] (Intercept) -1.174 0.378 -3.11 1.87e-03
## [ibex] AreaSB -1.257 0.664 -1.89 5.84e-02
## [ibex] AreaSU -0.606 0.552 -1.10 2.73e-01
## [ibex] AreaTB -0.933 0.648 -1.44 1.50e-01
## [ibex] EL 0.265 0.187 1.42 1.56e-01
## [ibex] SL 0.374 0.169 2.21 2.71e-02
## [leopardo] (Intercept) -1.407 0.427 -3.29 9.97e-04
## [leopardo] AreaSB 2.047 1.275 1.61 1.08e-01
## [leopardo] AreaSU 2.530 0.517 4.90 9.80e-07
## [leopardo] AreaTB -1.876 0.814 -2.30 2.12e-02
## [ibex:leopardo] (Intercept) 1.484 0.560 2.65 8.01e-03
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [ibex] (Intercept) -3.524 0.112 -31.47 2.53e-217
## [ibex] CT.sensmedium 0.686 0.227 3.02 2.51e-03
## [leopardo] (Intercept) -3.051 0.066 -46.26 0.00e+00
## [leopardo] CT.sensmedium -1.495 0.392 -3.81 1.39e-04
##
## AIC: 3918.218
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 104
## Bootstrap iterations: 04.5 Interazione leopardo - lupo
## model negLogLike nPars AIC delta AICwt
## 1 lup_leo_dip 1748.239 11 3518.478 0.00000 0.6472368
## 2 lup_leo_ind 1749.846 10 3519.692 1.21383 0.3527632In questo caso il modello più supportato contiene l’interazione. Modello senza interazione:
##
## Call:
## occuMulti(detformulas = c("~CT.sens", "~1"), stateformulas = c("~Area",
## "~EL+Dis", "~0"), data = dt5, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [leopardo] (Intercept) -0.725 0.303 -2.39 1.68e-02
## [leopardo] AreaSB 5.786 NaN NaN NaN
## [leopardo] AreaSU 2.263 0.464 4.88 1.07e-06
## [leopardo] AreaTB -2.166 0.789 -2.75 6.05e-03
## [lupo] (Intercept) -0.764 0.221 -3.46 5.35e-04
## [lupo] EL -0.299 0.205 -1.46 1.44e-01
## [lupo] Dis 0.506 0.198 2.55 1.07e-02
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [leopardo] (Intercept) -3.05 0.0661 -46.15 0.00e+00
## [leopardo] CT.sensmedium -1.83 0.2355 -7.79 6.67e-15
## [lupo] (Intercept) -3.95 0.1431 -27.61 8.13e-168
##
## AIC: 3519.692
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 62
## Bootstrap iterations: 0Modello con interazione:
##
## Call:
## occuMulti(detformulas = c("~CT.sens", "~1"), stateformulas = c("~Area",
## "~EL+Dis", "~1"), data = dt5, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [leopardo] (Intercept) -0.996 0.337 -2.95 3.16e-03
## [leopardo] AreaSB 1.375 1.020 1.35 1.78e-01
## [leopardo] AreaSU 2.265 0.465 4.87 1.10e-06
## [leopardo] AreaTB -2.075 0.791 -2.62 8.68e-03
## [lupo] (Intercept) -1.169 0.331 -3.53 4.15e-04
## [lupo] EL -0.231 0.205 -1.13 2.60e-01
## [lupo] Dis 0.408 0.205 1.99 4.62e-02
## [leopardo:lupo] (Intercept) 0.738 0.450 1.64 1.01e-01
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [leopardo] (Intercept) -3.05 0.066 -46.21 0.00e+00
## [leopardo] CT.sensmedium -1.45 0.390 -3.73 1.95e-04
## [lupo] (Intercept) -3.94 0.143 -27.62 6.00e-168
##
## AIC: 3518.478
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 46
## Bootstrap iterations: 04.6 Interazione stambecco - lupo
lup_ibex_dip <- occuMulti(detformulas = c('~1','~CT.sens'), stateformulas = c('~EL+Dis', '~Area+EL+SL', '~1'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))
lup_ibex_ind <- occuMulti(detformulas = c('~1','~CT.sens'), stateformulas = c('~EL+Dis', '~Area+EL+SL', '~0'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))## model negLogLike nPars AIC delta AICwt
## 2 lup_ibex_ind 1230.308 12 2484.616 0.0000000 0.5418256
## 1 lup_ibex_dip 1229.476 13 2484.951 0.3353889 0.4581744In questo caso il modello più supportato non contiene l’interazione. Modello senza interazione:
##
## Call:
## occuMulti(detformulas = c("~1", "~CT.sens"), stateformulas = c("~EL+Dis",
## "~Area+EL+SL", "~0"), data = dt1, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [lupo] (Intercept) -0.764 0.221 -3.462 0.000537
## [lupo] EL -0.299 0.205 -1.463 0.143355
## [lupo] Dis 0.506 0.198 2.552 0.010716
## [ibex] (Intercept) -0.682 0.307 -2.221 0.026351
## [ibex] AreaSB -0.576 0.490 -1.175 0.239802
## [ibex] AreaSU 0.174 0.424 0.411 0.681279
## [ibex] AreaTB -1.295 0.615 -2.106 0.035181
## [ibex] EL 0.270 0.187 1.448 0.147610
## [ibex] SL 0.368 0.169 2.172 0.029825
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [lupo] (Intercept) -3.950 0.143 -27.61 8.45e-168
## [ibex] (Intercept) -3.520 0.112 -31.55 1.94e-218
## [ibex] CT.sensmedium 0.684 0.227 3.02 2.54e-03
##
## AIC: 2484.616
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 107
## Bootstrap iterations: 0##
## Call:
## occuMulti(detformulas = c("~1", "~CT.sens"), stateformulas = c("~EL+Dis",
## "~Area+EL+SL", "~1"), data = dt1, control = list(maxit = 1000,
## trace = TRUE, REPORT = 1), maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [lupo] (Intercept) -0.599 0.258 -2.323 0.02017
## [lupo] EL -0.290 0.208 -1.391 0.16419
## [lupo] Dis 0.524 0.199 2.630 0.00855
## [ibex] (Intercept) -0.520 0.331 -1.572 0.11587
## [ibex] AreaSB -0.558 0.491 -1.135 0.25637
## [ibex] AreaSU 0.225 0.427 0.528 0.59766
## [ibex] AreaTB -1.329 0.614 -2.165 0.03035
## [ibex] EL 0.264 0.189 1.398 0.16209
## [ibex] SL 0.367 0.169 2.169 0.03011
## [lupo:ibex] (Intercept) -0.581 0.456 -1.274 0.20272
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [lupo] (Intercept) -3.955 0.144 -27.55 4.37e-167
## [ibex] (Intercept) -3.521 0.112 -31.53 3.46e-218
## [ibex] CT.sensmedium 0.683 0.227 3.01 2.61e-03
##
## AIC: 2484.951
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 98
## Bootstrap iterations: 0##Selezione dei modelli multispecie, a coppie, testando l’ipotesi di interazione
## model negLogLike nPars AIC delta AICwt
## 2 dom_ibex_ind 3319.739 17 6673.477 0.00000000 0.5093794
## 1 dom_ibex_dip 3318.776 18 6673.552 0.07504392 0.4906206
## 11 dom_leo_dip 3837.190 16 7706.380 0.00000000 0.7496418
## 21 dom_leo_ind 3839.286 15 7708.573 2.19340610 0.2503582
## 12 dom_lup_dip 3118.852 14 6265.704 0.00000000 0.7772587
## 22 dom_lup_ind 3121.102 13 6268.203 2.49952463 0.2227413
## 13 ib_leo_dip 1944.109 15 3918.218 0.00000000 0.9583079
## 23 ib_leo_ind 1948.244 14 3924.488 6.26971507 0.0416921
## 14 lup_leo_dip 1748.239 11 3518.478 0.00000000 0.6472368
## 24 lup_leo_ind 1749.846 10 3519.692 1.21383043 0.3527632
## 25 lup_ibex_ind 1230.308 12 2484.616 0.00000000 0.5418256
## 15 lup_ibex_dip 1229.476 13 2484.951 0.33538885 0.45817444.7 Riassunto grafico della selezione dei modelli per le interazioni inter-specifiche
Estimates for pairwise co-occurrence probabilities. Error bars show 95% confidence interval. Light blue indicates pairwise interactions supported by AIC, while orange indicates non supported pairwise interactions.
Estimates for pairwise co-occurrence probabilities. Error bars show 95% confidence interval. Light blue indicates pairwise interactions supported by AIC, while orange indicates non supported pairwise interactions.
4.7.1 predizioni multi-species
## $Livestock
## $Livestock$KS
## mean standard_dev SE
## 0.60534397 0.14784248 0.07206244
##
## $Livestock$SB
## mean standard_dev SE
## 0.44057281 0.15844494 0.08097191
##
## $Livestock$SU
## mean standard_dev SE
## 0.63971750 0.25437865 0.07258252
##
## $Livestock$TB
## mean standard_dev SE
## 0.29884588 0.12227894 0.07414067
##
##
## $Ibex
## $Ibex$KS
## mean standard_dev SE
## 0.34078750 0.10464736 0.08210319
##
## $Ibex$SB
## mean standard_dev SE
## 0.20373393 0.06154862 0.06897559
##
## $Ibex$SU
## mean standard_dev SE
## 0.37829737 0.11169770 0.08623092
##
## $Ibex$TB
## mean standard_dev SE
## 0.15006472 0.06132802 0.06992785
##
##
## $`Snow-leopard`
## $`Snow-leopard`$KS
## mean standard_dev SE
## 0.32015771 0.00000000 0.06640912
##
## $`Snow-leopard`$SB
## mean standard_dev SE
## 0.7382111 0.0000000 0.3053235
##
## $`Snow-leopard`$SU
## mean standard_dev SE
## 0.82314430 0.00000000 0.05122325
##
## $`Snow-leopard`$TB
## mean standard_dev SE
## 0.05237575 0.00000000 0.03617987
##
##
## $Wolf
## $Wolf$KS
## mean standard_dev SE
## 0.30501212 0.08089061 0.06422433
##
## $Wolf$SB
## mean standard_dev SE
## 0.3444104 0.1217007 0.0682746
##
## $Wolf$SU
## mean standard_dev SE
## 0.38605452 0.12440922 0.08421927
##
## $Wolf$TB
## mean standard_dev SE
## 0.25220791 0.06903757 0.06117556